Preparing the data

### Notes
# install.packages('ggplot2',dependencies = T)
# Test code line with command+enter
# update.packages(checkBuilt = T,ask = F)
# head(ds) # first few lines of a dataframe
# str(ds) #show structure of dataframe
###

# Setting variables:
options(width=109)
outT <- 2.5 # Outlier factor

# Loading libraries
library(ggplot2) # for plotting
library(plyr)

#Importing the data
#data_dir <- '/Volumes/PK/'
#homeDir <- path.expand('~/..')
homeDir <- path.expand('~')
dataDir <- paste(homeDir, 'Dropbox/Projects/fc/fc/data/usable/extracted', sep='/')
dataFiles <- dir(dataDir, pattern='csv')
nSubj <- length(dataFiles) # number of subjects

# Going through all subjects' data:
df <- data.frame() # data frame that combines all subjects' data into a single ds
#curSubjN <- 6 #temp

for(curSubjN in 1:nSubj){
    # Loading the data;
    ds <- read.csv(paste(dataDir, dataFiles[curSubjN], sep='/'))
    # Shorter names for variables:
    colnames(ds) <- c('subjId', 'domEyeR', 'threshStHi', 'threshStLo', 'thresh', 'trialN', 'sentId',
                      'sentPx', 'congr', 'fam', 'locTop', 'cued', 'corr', 'broken', 'st')
    # Non-blank data set:
    cds <- ds[ds$sentId<31,]
    # Removing the outliers:
    outlLo <- cds$st>(mean(cds$st)-outT*sd(cds$st)) # binary for every row
    outlHi <- cds$st<(mean(cds$st)+outT*sd(cds$st))
    cds$outlNlo <- sum(!outlLo) # counting the number of outliers before removing
    cds$outlNhi <- sum(!outlHi)
    cds <- cds[outlLo,] # removing the outliers
    cds <- cds[outlHi,]
    cds <- cds[!is.na(cds$subjId),] # usually don't need this; remove NA rows
    # Normalized suppression times:
    cds$stNorm <- cds$st / mean(cds$st)    
    # Binding the subject ds to the common data frame:
    df <- rbind(df, cds)
}
# Prettying up the data frame:
df$subjId <- factor(df$subjId)
df$Familiarity[df$fam==1] <- 'Familiar\n(Upright)'
df$Familiarity[df$fam==0] <- 'Unfamiliar\n(Inverted)'
df$Congruency[df$congr==1] <- 'Congruent'
df$Congruency[df$congr==0] <- 'Incongruent'
df$Attention[df$cued==1] <- 'Cued'
df$Attention[df$cued==0] <- 'Uncued'
# Binning
df$bin <- 6
df$bin[df$trialN<500] <- 5
df$bin[df$trialN<400] <- 4
df$bin[df$trialN<300] <- 3
df$bin[df$trialN<200] <- 2
df$bin[df$trialN<100] <- 1
head(df)
##   subjId domEyeR threshStHi threshStLo  thresh trialN sentId sentPx congr fam locTop cued corr broken     st
## 1     14       1    0.31181    0.23408 0.27294      1      5    133     1   0      1    1    1      1 1.4901
## 2     14       1    0.31181    0.23408 0.27294      2     28    108     1   0      0    0    1      1 1.2797
## 3     14       1    0.31181    0.23408 0.27294      5     27    104     0   0      1    1    1      1 1.5523
## 4     14       1    0.31181    0.23408 0.27294      6     21    109     0   0      1    1    1      1 1.4809
## 5     14       1    0.31181    0.23408 0.27294      8     18    117     0   0      1    1    1      1 1.3432
## 6     14       1    0.31181    0.23408 0.27294      9     10    128     0   0      0    0    1      1 1.1362
##   outlNlo outlNhi    stNorm            Familiarity  Congruency Attention bin
## 1       0       7 1.0650477 Unfamiliar\n(Inverted)   Congruent      Cued   1
## 2       0       7 0.9146645 Unfamiliar\n(Inverted)   Congruent    Uncued   1
## 3       0       7 1.1095051 Unfamiliar\n(Inverted) Incongruent      Cued   1
## 4       0       7 1.0584720 Unfamiliar\n(Inverted) Incongruent      Cued   1
## 5       0       7 0.9600510 Unfamiliar\n(Inverted) Incongruent      Cued   1
## 6       0       7 0.8120980 Unfamiliar\n(Inverted) Incongruent    Uncued   1

Quality control

Summary stats per participant

ddply(df, c('subjId'), summarise, outlNumLow = median(outlNlo), outlNumHigh = median(outlNhi),
      cuedCaught = sum(cued & corr), uncuedCaught = sum(!cued & corr))
##    subjId outlNumLow outlNumHigh cuedCaught uncuedCaught
## 1      14          0           7        116          116
## 2      17          0           9        115          116
## 3      19          0          10        113          107
## 4      21          0           6        118          116
## 5      23          0           6        117          116
## 6      24          0           0        105          101
## 7      25          0           9        115          116
## 8      26          0           9        118          113
## 9      27          0          11        110          115
## 10     28          1           7        115          116
## 11     29          0          14        114          112
## 12     30          0           5        118          117
## 13     31          0           7        116          115
## 14     32          0           7        119          114
## 15     33          0          11        107          104
## 16     34          0           9        113          117
## 17     36          0           6        118          115
## 18     37          0           4        119          117
## 19     38          1           4        115          116
## 20     39          0           7        117          116
## 21     41          0           7        119          114
## 22     42          0           6        120          114
## 23     43          4           4        115          116
## 24     44          0           6        118          111

sts across trials

for(curSubjN in 1:nSubj){
    curSubjId <- levels(df$subjId)[curSubjN]
    print(curSubjId)
    ss <- df[df$subjId==curSubjId,]
    p <- ggplot(data=ss, aes(x=trialN, y=st)) +  geom_line() + theme_bw() + ylim(.5,2) +
        labs(x='Trial Number', y='Suppression Time', title=paste('Participant', as.character(curSubjId))) + 
        theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(),
              axis.text=element_text(size=8), axis.title=element_text(size=9),
              legend.text=element_text(size=8), legend.title=element_text(size=9),
              legend.key = element_blank(), #legend.margin=unit(-.04, 'in'),
              legend.background = element_rect(fill='transparent'))
    plot(p)
}
## [1] "14"

## [1] "17"

## [1] "19"

## [1] "21"
## Warning: Removed 1 rows containing missing values (geom_path).

## [1] "23"
## Warning: Removed 1 rows containing missing values (geom_path).

## [1] "24"

## [1] "25"

## [1] "26"

## [1] "27"

## [1] "28"

## [1] "29"
## Warning: Removed 1 rows containing missing values (geom_path).

## [1] "30"

## [1] "31"

## [1] "32"

## [1] "33"

## [1] "34"

## [1] "36"

## [1] "37"

## [1] "38"

## [1] "39"
## Warning: Removed 2 rows containing missing values (geom_path).

## [1] "41"
## Warning: Removed 2 rows containing missing values (geom_path).

## [1] "42"

## [1] "43"
## Warning: Removed 1 rows containing missing values (geom_path).

## [1] "44"

Individual plots

Congruence X Familiarity X Cue

for(curSubjN in 1:nSubj){
    curSubjId <- levels(df$subjId)[curSubjN]
    ss <- df[df$subjId==curSubjId,]
    p <- ggplot(data=ss, aes(x=factor(Familiarity), y=st, colour=factor(Congruency))) +  geom_boxplot() + 
        facet_grid(~Attention) + theme_bw() +
        labs(x='Familiarity', y='Suppression Time', colour='Congruency',
             title=paste('Participant', as.character(curSubjId))) + 
        theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(),
              axis.text=element_text(size=8), axis.title=element_text(size=9),
              legend.text=element_text(size=8), legend.title=element_text(size=9),
              legend.key = element_blank(), #legend.margin=unit(-.04, 'in'),
              legend.background = element_rect(fill='transparent'))
    plot(p)
}

Group plots

Is cueing effective across participants?

p <- ggplot(data=df, aes(x=Attention, y=st)) + geom_boxplot() + facet_grid(~subjId) + theme_bw() +
    labs(x='Attention effect X Subject', y='Suppression Time', colour='Attention',
         title='Effect of attentional cueing') +
    theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(),
          axis.text=element_text(size=8), axis.title=element_text(size=9),
          legend.text=element_text(size=8), legend.title=element_text(size=9),
          legend.key = element_blank(), #legend.margin=unit(.0, 'in'),
          legend.background = element_rect(fill='transparent'))
plot(p)

Is cueing effective across trials?

p <- ggplot(data=df, aes(x=Attention, y=st)) + geom_boxplot() + facet_grid(~bin) + theme_bw() +
    labs(x='Attention effect X 100-trial bins', y='Suppression Time', colour='Attention',
         title='Effect of attentional cueing') +
    theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(),
          axis.text=element_text(size=8), axis.title=element_text(size=9),
          legend.text=element_text(size=8), legend.title=element_text(size=9),
          legend.key = element_blank(), #legend.margin=unit(.0, 'in'),
          legend.background = element_rect(fill='transparent'))
plot(p)

Congruence X Familiarity X Cue

dfSum <- ddply(df, .(subjId, Attention, Familiarity, Congruency), summarise, 
               `Mean Normalized ST` = mean(stNorm, na.rm=T))
p <- ggplot(data=dfSum, aes(x=factor(Familiarity), y=`Mean Normalized ST`, 
                            colour=factor(Congruency))) + 
    geom_boxplot() +  facet_grid(~Attention) + theme_bw() +
    labs(x='Familiarity', y='Suppression Time', colour='Congruency',
         title='Congruence X Familiarty X Cue') + 
    theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(),
          axis.text=element_text(size=8), axis.title=element_text(size=9),
          legend.text=element_text(size=8), legend.title=element_text(size=9),
          legend.key = element_blank(), #legend.margin=unit(-.04, 'in'),
          legend.background = element_rect(fill='transparent'))
plot(p)

p <- ggplot(data=dfSum, aes(x=factor(Familiarity), y=`Mean Normalized ST`, 
                            group=factor(subjId), colour=factor(subjId))) + 
    geom_line() + facet_grid(~Attention*Congruency) + theme_bw() +
    labs(x='Familiarity', y='Suppression Time', colour='Congruency',
         title='Congruence X Familiarty X Cue') + 
    theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(),
          axis.text=element_text(size=8), axis.title=element_text(size=9),
          legend.text=element_text(size=8), legend.title=element_text(size=9),
          legend.key = element_blank(), #legend.margin=unit(-.04, 'in'),
          legend.background = element_rect(fill='transparent'))
plot(p)

Statistical tests

Centered data set

df$congrC <- -1
df$congrC[df$congr==1] <- 1
df$famC <- -1
df$famC[df$fam==1] <- 1
df$cuedC <- -1
df$cuedC[df$cued==1] <- 1
df$topSubj <- 0
df$topSubj[df$locTop==df$cued] <- 1

Mixed linear regressions

library(lme4)
## Loading required package: Matrix
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
db <- '~/Dropbox/' # on Mac
source(paste(db, 'Prog/R/myFunctions/pvalfn.R', sep=''))
#mFull <- aov(st ~ fam * cued * congr + Error(subjId/(fam*cued*congr)), df)
#mFull <- aov(st ~ fam * cued * congr + Error(subjId), df)
#summary(mFull)
mFull <- lmer(st ~ famC * cuedC * congrC + (1|subjId/sentId), df)
pvalfn(mFull)
##                     estm    se       df   tVal  pVal star
## (Intercept)        1.506 0.038   23.013 39.273 0.000  ***
## famC              -0.001 0.006 5552.063 -0.167 0.867     
## cuedC             -0.037 0.006 5552.063 -6.116 0.000  ***
## congrC             0.003 0.006 5552.049  0.484 0.629     
## famC:cuedC         0.000 0.006 5552.052 -0.070 0.944     
## famC:congrC        0.004 0.006 5552.043  0.628 0.530     
## cuedC:congrC      -0.008 0.006 5552.047 -1.405 0.160     
## famC:cuedC:congrC -0.006 0.006 5552.051 -1.029 0.303

Bayes Factor

Full ‘top’ model results

Numbers above 10 indicate a significant improvement of the model (compared to the full model) when a term is /dropped/.

library(BayesFactor)
## Loading required package: coda
## ************
## Welcome to BayesFactor 0.9.12-2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
## 
## Type BFManual() to open the manual.
## ************
df$subjId <- as.factor(df$subjId)
df$sentId <- as.factor(df$sentId)
df$fam<- as.factor(df$fam)
df$cued<- as.factor(df$cued)
df$congr<- as.factor(df$congr)
bfTop <- anovaBF(st ~ fam * cued * congr + subjId + sentId, data=df, 
                 whichRandom = c('subjId','sentId'), whichModels='top')
bfTop
## Bayes factor top-down analysis
## --------------
## When effect is omitted from congr + fam + cued + congr:fam + congr:cued + fam:cued + congr:fam:cued + subjId + sentId , BF is...
## [1] Omit congr:cued:fam : 13.38371     <U+00B1>19.05%
## [2] Omit cued:fam       : 24.07596     <U+00B1>4.76%
## [3] Omit congr:cued     : 10.55362     <U+00B1>10.21%
## [4] Omit congr:fam      : 21.29736     <U+00B1>5.57%
## [5] Omit cued           : 2.547776e-07 <U+00B1>5.03%
## [6] Omit fam            : 35.13981     <U+00B1>5.47%
## [7] Omit congr          : 33.30138     <U+00B1>8.11%
## 
## Against denominator:
##   st ~ congr + fam + cued + congr:fam + congr:cued + fam:cued + congr:fam:cued + subjId + sentId 
## ---
## Bayes factor type: BFlinearModel, JZS

Improvement over minimal cued-only model?

Numbers below 1 indicate a significant detriment to the cued-only model if a given term is included.

bfFull<- anovaBF(st ~ fam * cued * congr + subjId + sentId, data=df, 
                 whichRandom = c('subjId','sentId'))
bfCued <- anovaBF(st ~ cued + subjId + sentId, data=df,
                  whichRandom = c('subjId','sentId'))
bfFull/ bfCued
## Bayes factor analysis
## --------------
## [1] congr + subjId + sentId                                                                    : 8.91964e-09  <U+00B1>2.01%
## [2] fam + subjId + sentId                                                                      : 8.196144e-09 <U+00B1>1.64%
## [3] congr + fam + subjId + sentId                                                              : 2.595342e-10 <U+00B1>1.77%
## [4] congr + fam + congr:fam + subjId + sentId                                                  : 1.344002e-11 <U+00B1>2.64%
## [5] cued + subjId + sentId                                                                     : 1            <U+00B1>0%
## [6] congr + cued + subjId + sentId                                                             : 0.03391136   <U+00B1>3.2%
## [7] fam + cued + subjId + sentId                                                               : 0.03079315   <U+00B1>2.55%
## [8] congr + fam + cued + subjId + sentId                                                       : 0.0009933776 <U+00B1>3.31%
## [9] congr + fam + congr:fam + cued + subjId + sentId                                           : 5.280046e-05 <U+00B1>5.87%
## [10] congr + cued + congr:cued + subjId + sentId                                               : 0.004026464  <U+00B1>4.84%
## [11] congr + fam + cued + congr:cued + subjId + sentId                                         : 0.0001077874 <U+00B1>2.61%
## [12] congr + fam + congr:fam + cued + congr:cued + subjId + sentId                             : 6.121374e-06 <U+00B1>10.36%
## [13] fam + cued + fam:cued + subjId + sentId                                                   : 0.001219005  <U+00B1>2.35%
## [14] congr + fam + cued + fam:cued + subjId + sentId                                           : 4.12929e-05  <U+00B1>2.72%
## [15] congr + fam + congr:fam + cued + fam:cued + subjId + sentId                               : 2.104796e-06 <U+00B1>3.72%
## [16] congr + fam + cued + congr:cued + fam:cued + subjId + sentId                              : 5.206666e-06 <U+00B1>6.31%
## [17] congr + fam + congr:fam + cued + congr:cued + fam:cued + subjId + sentId                  : 2.402968e-07 <U+00B1>5.59%
## [18] congr + fam + congr:fam + cued + congr:cued + fam:cued + congr:fam:cued + subjId + sentId : 2.339387e-08 <U+00B1>5.77%
## 
## Against denominator:
##   st ~ cued + subjId + sentId 
## ---
## Bayes factor type: BFlinearModel, JZS

Bayesian Regression

regFull <- regressionBF(st ~ famC * cuedC * congrC, data=df)
regFull
## Bayes factor analysis
## --------------
## [1] famC                  : 0.03084727 <U+00B1>0%
## [2] cuedC                 : 206713.5   <U+00B1>0%
## [3] congrC                : 0.03252247 <U+00B1>0%
## [4] famC + cuedC          : 9760.311   <U+00B1>0%
## [5] famC + congrC         : 0.00157279 <U+00B1>0.01%
## [6] cuedC + congrC        : 10366.37   <U+00B1>0%
## [7] famC + cuedC + congrC : 622.0742   <U+00B1>0.01%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS
regCued <- regressionBF(st ~ cuedC, data=df)
regFull / regCued
## Bayes factor analysis
## --------------
## [1] famC                  : 1.492272e-07 <U+00B1>0%
## [2] cuedC                 : 1            <U+00B1>0%
## [3] congrC                : 1.573311e-07 <U+00B1>0%
## [4] famC + cuedC          : 0.04721661   <U+00B1>0%
## [5] famC + congrC         : 7.608549e-09 <U+00B1>0.01%
## [6] cuedC + congrC        : 0.05014849   <U+00B1>0%
## [7] famC + cuedC + congrC : 0.003009355  <U+00B1>0.01%
## 
## Against denominator:
##   st ~ cuedC 
## ---
## Bayes factor type: BFlinearModel, JZS